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FORTRAN IV SUBROUTINES FOR COUPLING COEFFICIENTS AND 


MATRIX ELEMENTS IN THE QUANTUM MECHANICAL 
THEORY OF ANGULAR MOMENTUM 
by William F. Ford and Richard C. Braley 
Lewis Research Center 

SUMMARY 

Subroutines written in FORTRAN IV are presented for calculating Clebsch-Gordan 
coefficients, Racah coefficients, 9-j coefficients, reduced rotation matrix elements, and 
other related quantities. Considerable attention is paid to matters of speed and accura- 
cy, and the coding is designed to resemble the corresponding mathematical equations as 
closely as possible. 


INTRODUCTION 

The investigation of quantum mechanical systems of two or more particles generally 
involves the coupling of the individual angular momenta to some total angular momentum. 
The study of this process leads quite naturally to the Clebsch-Gordan coefficients, which 
describe the coupling of a pair of particles, and to the Racah and other coefficients, 
which describe the transformation from one coupling scheme to another. 

Extensive tabulations of such quantities have been available for many years, and 
these are very useful for calculations which one might undertake with the aid of a desk 
calculator. With the advent of high-speed electronic computers, however, the use of ex- 
tensive tabulations becomes more of a hindrance than a help, and one is faced with the 
need to generate the relevant coefficients as and when they are required. This is partic- 
ularly true of the reduced matrix elements of the rotation operator, which occur in many 
applications, because they depend on a continuous variable. 

A general-purpose routine designed to furnish some quantity in widespread use, say 
for example, a Clebsch-Gordan coefficient, ought to satisfy certain requirements. It 
should be accurate, of course, and fast; in addition it must be easy to use and designed 
in such a way that the user can modify it readily to suit a particular application. 


A common problem arising in the evaluation of the Ciebsch-Gordan and Racah coef- 
ficients and rotation matrix elements is the appearance of large factorials in the numer- 
ators and denominators of the terms in the series. Direct use of the factorials may lead 
to roundoff or overflow, while attempts to avoid the problem by introducing binomial co- 
efficients are clumsy and time-consuming. The method' used here involves logarithms 
of the various factorials, which are precalculated and stored in an array; the exponent- 
iation is not performed until cancellation occurs via subtraction of the logarithms. 

In the following sections we present techniques for computing the Ciebsch-Gordan 
coefficient, the Racah coefficient, the reduced rotation matrix element, and the 9-j co- 
efficient. FORTRAN IV subroutines which employ these techniques are given in the ap- 
pendices. The subroutine for evaluating the 9-j coefficient also contains an entry for 
computing the reduced matrix element of the spin-angle tensor product. 

Within each subroutine the coding has been designed to resemble the actual equa- 
tions, so that a potential user can easily make modifications if desired. In order to 
maintain this close correspondence and at the same time produce efficient coding, the 
basic equations have been rearranged slightly from the forms given by Brink and Satchler 
(ref. 1). The reader should also note that the same symbol may be used to represent 
different quantities in different sections, and so the definition given in one section ap- 
plies to that section only. 


CLEBSCH-GORDAN COEFFICIENTS 

The expression given by Brink and Satchler (ref. 1, p. 34) for the Ciebsch-Gordan 
coefficient may be written 


<il m l, j 2 m 2| JM > =CG 


where 


C = A(j 1 j 2 J) [(jj - m^Kjg +m 2 )’.(j 1 + m 1 )’.(j 2 - m 2 )’.(J +M)I(J - M)l] ^ 2 

and 

G = |/2J + l»l) n [(jj +j 2 - J - n)’. (jj - mj - n)».(j 2 +m 2 - n)l (J - j 2 +mj +n)l 
n 

x (J - j 1 - mj +n)ln’.] -1 
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Here the symbol A, which will appear in the next section as well, stands for a function 
defined by 

A(abc) = [(a +b - c)!(b + c - a)l (c + a - b)'.] 1 / 2 [(a +b + c + l)'.]*/ 2 

The appearance of factorials in the numerators and denominators of the above ex- 
pressions makes it undesirable to evaluate them as written. Our philosophy here and in 
the remaining sections will be to avoid the direct use of factorials by means of loga- 
rithms, and to arrange matters so that the first term in the series has the value unity. 
Thus, to evaluate C, we make use of the fact that 

nl = T(n + 1) = exp {log T(n + 1)} 

to write 

9 

C = exp^i P = XCNj) - X(N 1q ) 

where X(N) s log T(N) and 

Nj = 1 + (jj + j 2 - J) N g = 1 + (jj + mj) 

N 2 = 1 + (j 2 + J - J x ) N ? = 1 + (j 2 - m 2 ) 

N 3 = 1 +(J + jj - j 2 ) N g = 1 + (J + M) 

N 4 = l+(j 1 -m 1 ) N g = 1 + (J - M) 

Ng = 1 + (j 2 + m 2 ) Njq = 1 + (jj + ^2 + J + 1) 

= Nj +N 2 +Nj - 1 

The are always integers, and so C vanishes if any is smaller than unity; this 
condition is equivalent to all the usual triangle inequalities and also the z- component 
inequalities for the Clebsch-Gordan coefficient. Because X(N) is needed only for integer 
values of N, it can be conveniently precalculated and stored as an array labeled X. 

To evaluate G, we write 

n 2 

G = V G(n) = G(n x ) £ H(n) 
n n=n^ 
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where H(n^) = 1, 


H(n) = H(n - 1) - G ^. 

G(n - 1) 

and n-^ (n 2 ) is the minimum (maximum) value attained by n. The recursion relation for 
H(n) reduces to 


(n - Kj)(n - K 2 )(n - Kg) 



(n - K 4 )(n - Kg)n 

where 


K i = N i 

K 4 = N 4 - N 3 

k 2 = n 4 

K 5 = N 5 “ N 2 

K 3 = N 5 


It is easy to verify that n 2 is the smallest of (K 4 ,K 2 ,Kg) less one, and n^ the largest 
of(0,K 4 ,K 5 ). 

The first term G(n.), which has been factored out of the series, can be evaluated 

using the same technique as for C: 


n . 

6 

G(n x ) = j/2J + 1 (-1) 1 exp(-Q) Q = £ X ( K J) 

i=l 

where 


K i = K i - n i 

K 4 = n i + 1 " K 4 

K 2 = K 2 ' n l 

K^nj + 1-K 5 

k 3 = k 3 - n 1 

K e = n l + 1 


The FORTRAN IV subroutine CG which evaluates the Clebsch-Gordan coefficient by 
means of these equations is presented in appendix A. The coding has been arranged so 
as to closely resemble the actual equations. The subroutine is written as a function 
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subprogram, so that a typical call might be 


SUM = SUM + ALPHA * CG(J1, Ml, J2, M2, J, M) 

The arguments Jl, . . . , M in the calling vector are to be supplied as floating point 
variables, so that either integer or half-integer values may be easily handled. 

Provision has also been made for evaluation of the Wigner 3-j coefficient, which is 
related to the Clebsch-Gordan coefficient by 


' j l h J \_(-l) 


V m l m 2 M j 


h-h 


io-M 


|/2J + 1 

A typical call for this function might be 


<j 1 m 1 J 2 m 2 J “ M ) 


SUM = SUM + ALPHA * THREEJ(J1, J2, J, Ml, M2, M) 


RACAH COEFFICIENTS 

The expression given by Brink and Satchler (ref. 1, p. 43) for the Racah coefficient 
may be written 


W(abcd; ef) = CG 


where 


C = A(abe) A(acf) A(bdf) A(cde) 

and 

G = (~ !) n ( a + b + c + d + 1 - n)I [(a + b - e - n)l (a + c - f - n)'. (b + d - f - n)l 

n 

x (c + d - e - n)I (e + f - a - d + n)i (e + f - b - c + n) , .n'.]”' 1 ‘ 

As in the preceding section, C can be evaluated more conveniently and accurately in the 
form 
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c 


= exp 



where X(N) = log r(N) and 


4 

P = ^ [ X < J i) +X < K i) + X ^ L i) - x ( M i)] 


i=l 


Jj = 1 + (a + b - e) 

K-^ = 1 + (b + e - a) 

Lj = 1 + (e + a - b) 

Mj = Jj +K X + Lj - 1 


Jg = 1 + (a + c - f) 
K 2 = 1 + (c + f - a) 
L 2 = 1 + (f + a - c) 

m 2 = J 2 + k 2 + L 2 
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J3 = 1 + (b + d - f) 

K 3 = 1 + (d + f - b) 

L 3 = 1 + (f + b - d) 

m 3 = J 3 + K 3 + L 3 - 1 


J 4 = 1 + (c + d - e) 

= 1 + (d + e - c) 

= 1 + (e + c - d) 

m 4 = J 4 + k 4 + l 4 - 1 


All of the quantities J^, K^, L-, and are integers, and the usual triangle inequal- 
ities on the arguments of the Racah coefficient are equivalent to requiring that all these 
integers (except for Mj) be greater than zero. 

To evaluate G we proceed as for the Clebsch-Gordan coefficient, writing 

n 2 

G "Yj = G ^l^ Yj 

n n=n^ 

with = 1 as before. In this case the recursion relation for H(n) reduces to 


(n - Ji)(n - J 2 )(n - Jo)(n - J 4 ) 

H(n) = H(n - 1) ± * 4 

(n - J 5 )(n - J g )(n - J 7 )n 
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where 


J 5 ~ J 1 ' L 2 


J 6 - j i - Lg 


Jrj — “ 1 


The first term G(nj) can be written 

G(n 1 )=(-1) ^xpC-Q) 


Q = Y, X(J i) “ X(J 8> 

i=l 


where 


J 1 = J 1 - 

n l 

J 5 - n 

f 

J 2 = J 2 " 

n l 

J 6 = n 

CO 

II 

CO 

•“3 

n l 

» 

Jrj = n 

J 4 “ J 4 ' 

n l 

= J, 


5 

*6 


The FORTRAN IV subroutine RACAH which evaluates the Racah coefficient by means 
of these equations is presented in appendix B. The coding has been arranged so as to 
closely resemble the actual equations. The subroutine is written as a function subpro- 
gram, so that a typical call might be 


SUM = SUM + ALPHA * RAC AH (A, B, C, D, E, F) 


The arguments A, . . . , F in the calling vector are to be supplied as floating point var- 
iables, so that either integer or half- integer values may be easily handled. 

Provision has also been made for evaluation of the normalized Racah coefficient, or 
U- coefficient, which is related to the Racah coefficient by 


U(abcd; ef) = y/(2e + l)(2f +1) W(abcd; ef) 

The U-coefficient has the convenient property that U = 1 whenever a, b, c, or d is 
zero, provided the triangular inequalities are satisfied. A typical call for this function 
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might be 


SUM = SUM + ALPHA * UCOEF(A, B, C, D, E, F) 

There is also an entry for evaluation of the Wigner 6-j coefficient, which is some- 
times used in preference to the Racah coefficient. .It is related to the Racah coefficient 
by 


/ abc l = (-l) a ’*' Ki+e w(abed; cf) 
IdefJ 


A typical call for this function might be 


SUM = SUM + ALPHA *SIXJ(A, B, C, D, E, F) 


REDUCED ROTATION MATRIX ELEMENTS 

The expression given by Brink and Satchler (ref» 1, p. 22) for the elements of the 
reduced rotation matrix may be written 


m-jing^) = CG 


where 


C = [(j + m j) 1 ( j - m 2 )'.(j - m^lQ + m 2 )'.] 1//2 


G ~ cos 




2n+m«-m 1 

1 [(j +mj - n)\ (j - m 2 - n)! 


» n*l - 1 


x (m 2 - m^ + n)’n’.J 


As in the preceding sections, C can be evaluated more conveniently and accurately in 
the form 


C = exp - P 
\ 2 . 


£ x( N[ ) 

i=l 
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where X(N) = log r(N) and 


Nj = 1 + (j + N 3 = 1 + (j - mj) 

N 2 = 1 + (j - m 2 ) N 4 = 1 + (j + m 2 ) 

As before, C vanishes if any of the integers is less than unity. 

Before evaluating G we must discuss its dependence on the rotation angle 9. Al- 
though G is periodic in 9 with a period of for half-integral j and a period of 2ir 
for integral j, in most applications the range of 9 is limited to 0 ^ 9 — it. We will 
take this to be the standard case, and discuss the exceptions later. 

Because the expression for G is a power series in tan — 9 , it is desirable to ar- 

1 i 2 

range matters so that 9 ^ — it. When 9 is larger than —n, this can be accomplished 

2 2 

by means of the relation 


dj 

m l m 2 


( e ) 


j- m i 

(- 1 ' Vm 2 <” ' e) 


We then write 


n 2 

G = G(n^) ^ H ( n ) 
n=n^ 

with H(nj) = 1 as before. In this case the recursion relation for H(n) reduces to 


. \2 (n - K. )(n - K 9 ) 

H(n) = - H(n - 1) tan - 6») 1 — 

2 / (n - K 3 )n 


where = Nj, K 2 = N 2 , and K 3 = - N^. 

The first term G(n^) can be written 


G( n i) =* ^cos i e'j ^tani-0^ (-1) 1 exp(-Q) 


Q = £ x(n’) 

i=l 


9 


where N = 2n^ - (nij - mg) and 


N l= K l- n l 

N 3 

= nj + 1 - 

» 

? 


N 2 = K 2 “ n l 

N 4 

= nj + 1 


Finally, let us consider the nonstandard case, where 9 is negative or larger than 
77 . For negative angles we may use the relation 


di 

m l m 2 


(e) = (-i) 


m r m 2 d j 


m l m 2 


(-9) 


For positive angles we may immediately reduce 9 to its value modulo 477. If 9 is 
then larger than 277 , we may use 


d j 

m l m 2 


(9) = - 2 *> 


and if 9 is then larger than tt, we may use 


m i m 2 


(0) = (-l) 


2j+mj-m 


2 dj 

m l m 2 


(277 - 9) 


The FORTRAN IV subroutine DJMM which evaluates the reduced rotation matrix 
element by means of these equations is presented in appendix C. The coding has been 
arranged so as to closely resemble the actual equations. The subroutine is written as a 
function subprogram, so that a typical call might be 

SUM = SUM + ALPHA *DJMM(M1, M2, J, THETA) 

The arguments Ml, . . . , THETA in the calling vector are to be supplied as floating 
point variables, so that either integer or half-integer values of the angular momentum 
quantum numbers may be easily handled. 
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9-j COEFFICIENTS AND RELATED QUANTITIES 


The expression given by Brink and Satchler (ref. 1, p. 46) for evaluating the Wigner 
9-j coefficient in terms of Racah coefficients may be written 


c i i b i* r 
4 l 2 s 2h 

Ll S J j 


hZ (2k + 1)W(Z 1 I 2 JS; LkjWtigZgSjS; s 2 k)W(Z jS^; jjk) 
k 


The summation is restricted by triangular inequalities such that the quantity 


A (l j Jk) A(z 2 Sk) A(j 2 s jk) 

does not vanish. By specializing one or another of the arguments it is possible to derive 
simpler expressions for the 9-j coefficient, but in general this will not be attempted 
here. 

The case where s-^ = s 2 = 1/2 occurs so frequently, however, that it is deserving 
of special attention. Only two values of S are possible here, namely S = 0 and S = 1. 
When S = 0 the 9-j coefficient may be obtained directly from the above series, since 
only one term survives; the result may be simplified using 


W(abcO; ef) = 


6 bf 5 ce 

f/(2e + l)(2f + 1) 


The general relation 


^T(2S + ljWfcjSjLJ; Sk) 
s 



/=W(j 1 s 1 Z 2 L; 1 1 k)W(Z 2 s 2 j 1 J; jgk) 


may then be applied to obtain the 9-j coefficient with S = 1 in terms of the 9-j coef- 
ficient with S = 0. The final result may be put in the compact form 


r i . > 

l iih 


Q 


LSJ 




Ll s j J 


> = 


A5 lj ~ B6 si 
C5 so _ D6 si 


li 



where 


A=w^iI J j 2 ; 


B = (4J +2)WU 2 ±j 1 J; j 2 EJw/j 1 ll 2 L; l ^ 


C = |/2(2J + 1) 


1 1 


D = 3(4J + 2)Wl LJ; lEj 


E = — (L + J +6 j^j) 
2 


The 9-j coefficients play an important part in evaluating the matrix elements of 
products of tensor operators 0 A case of special interest involves the spin-angle tensors 
defined by 


< m L* SM sl JM > Y L L °r S 


M L M S 


Here Y^ is a spherical harmonic and o^ 1 is a rank-S tensor in the spin- space of a 
particle with intrinsic spin. For spin- 1/2 particles ct q is the unit operator and aj is 
twice the spin operator, and the reduced matrix element turns out to be expressible as 

< l l\ hI^LSjll Z 2^2> s ^ <j 2 m 2’ JM lh m l>< Z l^H m ll^LSjl Z 2^2 m 2> 

m l m 2 

= |/2(2z 1 + l)(2j 2 + l)(2S + 1)(2J + 1)<Z i|| Y L |Z 2 )Q LSJ 


Here a further simplification is possible because of the appearance of the reduced matrix 
element of Y L , which is known to vanish unless l ^ + l 2 + L is even. After some alge- 
bra the result may be written 

^LSjIhf S 2 > = V 1 ? I< ' 1)J(il 2’ 
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where 


G LOL ~ 1 


'L1L 


X 1 ~ X 2 

/j(J + 1) 


_ X 1 +X 2 +A 
G L1L+1 " 

yA(2J + 1) 


with A = (L + J + l)/2 and X = (l - j)(2j + 1). 

The FORTRAN IV subroutine NINEJ which evaluates the 9-j coefficient in the gen- 
eral case is presented in appendix D. The coding has been arranged so as to closely re- 
semble the actual equations. The subroutine is written as a function subprogram, so 
that a typical call might be 


SUM = SUM + ALPHA * NINEJ(L1, SI, Jl, L2, S2, J2, L, S, J) 


The arguments LI, . . . , J in the calling vector are to be supplied as floating point 
variables, so that either integer or half-integer values may be easily handled. 

An entry has been provided for obtaining the 9-j coefficient in the special case 
where = Sg = 1/2. A typical call for this function might be 

SUM = SUM + ALPHA *QLSJ(L1, Jl, L2, J2, L, S, J) 

An entry has also been provided for obtaining the reduced matrix elements of the spin- 
angle tensor T . A typical call for this function might be 

SUM = SUM + ALPHA * TLSJ(L1, Jl, L2, J2, L, S, J) 


CONCLUDING REMARKS 

Techniques for computing the Clebsch-Gordan coefficient, the Racah coefficient, the 
reduced rotation matrix element, and the 9-j coefficient have been presented. The dif- 
ficulties associated with large factorials in the numerators and denominators of terms in 
the series were avoided by introducing the logarithms and by arranging matters so that 
the first term in the series had the value unity. Various related quantities were also de- 
fined and methods given for their computation. 
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FORTRAN IV subroutines based on these techniques are presented in the appendices. 
The coding has been designed to resemble the equations in the text, and although this 
leads to a slight loss in efficiency, it makes it easy for the user to make changes to suit 
his particular application. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, November 4, 1970, 

129-02. 
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APPENDIX A 
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SUBROUTINE FOR EVALUATING CLEBSCH-GORDON COEFFICIENT 


FUNCTION CG(J1*M1*J2?M2»J»M) 

LOGICAL FIRST/. TRUE./, ODD 
REAL J1 , Ml , J2,M2, J,M 
DOUBLE PRECISION SUM 
COMMON /CGCOM/X ( 200 ) 

DATA ZER0/5.0E-6/ 

0D0(N)=M0Q(N, 2) .EQ. 1 

C THIS SUBROUTINE CALCULATES THE C LE BSCH-GORDAN COEFFICIENT <J 1 , M 1, J 2, M 2/ J , M> 

MODE=l 
P = M 

GO TO 1 

ENTRY THREEJI J1 , J2, J, Ml, M2,M) 

M0DE=2 
P - — M 

1 IF (FIRST) GO TO 6 
C BEGIN CALCULATION 

2 CG=0. 0 

N = 0 . 1 +A B S ( MI+M2-P ) 

IF (N.NE.O) RETURN 

Nl=l.l+( J1+J2-J) 

N2 = l • 1 + ( J 2 + J— J 1 ) 

N3=1.1M J+JI-J2) 

N4=l.l+< Jl-Ml) 

N5=l.l+( J2+M2) 

N6=l. 1M Jl+Ml) 

N7- 1 . l + ( J 2-M 2 ) 

NR = 1.1 + (J + M) 

N9= l • 1+ ( J— M ) 

IF ( M I NO ( N l , N2 , N3 , N 4, N 5 , N6, N7» N8 , N9 ) .LT. 1 ) RETURN 
FAC TGR= 1 • 0 

IF (MODE.EQ.2.AND.ODD{Nl-N8) ) FACTOR =- FACTOR 

C TEST FOR SPECIAL VALUES 

IF ( Jl.LE. ZERO. OR. J2.LE. ZERO) GO TO 5 

IF ( Nl.EQ. 1 . AND. ( N8.EQ.1.0R.N9 .EO.l) ) GO TO 5 

C CONTINUE CALCULATION 

IF (MODE. EO.l) FACTO R=SQRT(2.0*J+1.0 ) *F ACTOR 
N=N1+N2+N3-1 

P-X(N1 ) + X ( N2 ) +X(N3 )+X(N4)+X(N5)+X(N6)+X(N7)+X (N8)+X(N9)-X(N ) 

K 1 = N 1 
K2 = N4 
K3 = N5 
K4=N4— N3 
K 5 = N5— N 2 

NMAX=MIN0(Kl f K2,K3)-l 
NMIN-MAXO (K4,K5,0 ) 

NP1=NMI N+ 1 

IF ( ODD ( NM I N ) ) FACTORY-FACTOR 
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5800 
5900 
6000 
6100 
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92 00 


TE RM=F ACTOR 
SUM=TERM 

IF (NMIN.EQ.NMAX1 GO TO 4 

DO 3 N=NP1 v NMAX 
0NE=(N-K1 )*(N-K2)*(N-K3) 

TW0=(N-K4)*IN-K5)*N 
TE RM=TE RM*ONE/TWO 

3 SUM=$UM+TFRM 

4 N1 =K1 — NMI N 
N2=K 2-NM IN 
N3 = K3 — NMI N 
N4=NP1-K4 
N5=NP 1— K5 

0=X (N1)+X(N2)+X(N3)+X(N4)+X(N5)+X(NP1) 

CG=EXP( 0 * 5*P— 0 ) *SUM 

IF (ABS(CG) .LE.ZEROI CG=0*0 

RETURN 

C SPECIAL CASES 

5 IF (M0DE.EQ.2) FACTOR=FACTOR/SQRT< 2* 0* J+ 1. 0) 
CG=F ACTOR 

RETURN 

C ARRAY X(N) = LOG(GAMMA(NM 

6 F I RST=. FALSE* 

DO 7 N= 1 » 200 
Q=N 

7 X( N)=ALGAMA ( Q) 

GO TO 2 

END 
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APPENDIX B 


SUBROUTINE FOR EVALUATING RACAH COEFFICIENT 
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1000 
lino 
1200 
13-00 
1400 
15 00 
1600 
17 00 
1800 
1900 
2000 
2100 
2200 
2300 
2400 
2 5 n 0 
2600 
2700 
2800 
2900 
3000 
3100 
3200 
3300 
3400 
3500 
3600 
3700 
3800 
3900 
40OQ 
4100 
4200 
4300 
4400 
4500 
4600 
4 7C0 
4800 
4900 
5000 
5100 
5200 
5300 
5400 
5500 
5600 
5700 
5800 


FUNCTION RACAH< A,B f C,D,E, F ) 

LOGICAL FIRST/. TRUE. /, 000 
DOUBLE PRECISION SUM 
COMMON /CGCOM/ X ( 200 ) 

DATA ZERO/5. OF-6/ 

ODD ( N ) = MOD (Nt2).FQ.l 

C THIS SUBROUTINE CALCULATES THE RACAH COEFFICIENT W ( A , B ,C f D ; E , F ) 

MODE = 1 
GO TO 1 

ENTRY UCOEF(A T B ? C f D,E,F) 

MODE = 2 
GO TO l 

F N TRY SIX J( AtB,F,D,C, F ) 

MODE =3 

1 IF (FIRST) GO TO 6 
C BEGIN CALCULATION 

2 R A C AH=0 • 0 
Jl=l. 1+ ( A* R-E ) 

K1-1.1 + (B+E-A ) 

L 1 = 1 . l+(E+4~B) 

IF (MINO( JlyKlvLl) .LT.l ) RETURN 

J2 = 1 . 1 + ( A* C-F ) 

K 2= 1 . 1MC + F-A) 

L?= 1 • 1 + ( F + A-C ) 

IF (MIN0(J?tK2tL2) .LT.l ) RETURN 

J3= 1 • 1+ ( B+D-F ) 

K. 3-l.l + (D + F- B) 

L 3= 1 . 1+ ( F + B-D ) 

IF ( MINO ( J3, K3, L3 ) .LT . 1 ) RETURN 

J 4= 1 • 1 + ( C + D— F ) 

K 4 - 1 .1 + 1 D+ E-C ) 

L4= 1 . 1 + ( E+C-D ) 

IF (MINCH J4,K4,L4) .LT.l) RETURN 

FAC T0R = 1 • 0 
M4=J4+K4+L4~ 1 
J7 - J 1 +M4— 1 

IF ( MODE. E 0. 3. ANO.ODDC J7) ) F AC TOR= — FACTOR 
C TEST FOR SPECIAL VALUES 

IF ( A. GT.ZFRO. AND. B.GT. ZERO. AND.C.GT. ZERO. AND. D.GT. ZERO) GO TO 3 
IF ( MODE. NE • 2 ) FACT OR = 0.5*F AC TOR/S QRT( ( E + 0 . 5 ) *1 F+0 • 5 M 
R AC AH = F AC TOR 
RETURN 

C CONTINUF CALCULATION 

3 P=X( Jl) + X( J2)fX( J3) + X( J4) *X(Klf +X(K2)+X (K3)+X( K4)+XtLi )+X(L2 J+X (L3I + XIL4) 
Ml= Jl+Kl+Ll-l 

M2=J2+K2+L2-1 
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5900 
6000 
6100 
6200 
6300 
6600 
6500 
6600 
6700 
6 8 00 
6900 
7000 
7100 
7200 
7300 
7600 
7500 
7600 
7700 
7800 
790 ^ 
8000 
8100 
8200 
8300 
8600 
8500 
8600 
8700 
8800 
8900 
9000 
91 00 
9200 
9300 
9600 
9500 
9600 
9700 
9800 


M3=J3+K3+L3-1 

P^P-X (Ml ) - X ( M2 ) — X ( M3 )— X ( M6 ) 

J 5 = J 1-L 2 
J6 = J 1 -L 3 

NM A X = M IN0(JlyJ2fJ3« J6 )-l 
NMIN=MAX0( J5t J6,0) 

NP 1 = NMT N + l 

IF (OQD(NMIN)) FACTORY-FACTOR 

TERM=FACTHR 

SU M=T F RM 

IF (NMIN.EQ.NMAX) GO TO 5 
DO 6 N=NP1,NMAX 

ONF = ( N — J 1 )*< N-J2 ) * { N— J 3)*(N— J6) 
TWD=(N-J5)*(N-J6)*(N-J7)*N 
TFPMyTERM«ONF/TWO 
6 SUM=SUM+TERM 

5 J 1 = J 1 — NM I N 
J2 = J? — NMI N 
J 3 = J3-NMT N 
J4= J6-NMIN 
J5 = NP 1 - J 5 
J6 = NP 1- J 6 
J7= J7-NMIM 

0=X( Jl)+X( J2 > + X( J3 )+XU6)+X ( J5)+X( J6)+X( NP1)-X(J7) 

RACAH=EXP(0.5*P-Q)*SUM 

IF ( ARS(RACAH) .LF^ZERO) RACAH=0.0 

RETURN 

C ARRAY X( N)=LOG(GAMMA(N) ) 

6 FIRST^FALSE. 

DO 7 N = I » 200 
Q=N 

7 X (N)= ALGAMA(G) 

GO tq 7 

END 



APPENDIX C 


SUBROUTINE FOR EVALUATING REDUCED ROTATION MATRIX ELEMENT 


100 

200 

300 

400 

500 

600 

700 

aoo 

900 
1 000 
1100 
1200 
1300 
1400 
1500 
1600 
1700 
1800 
1900 
2000 
2100 
2200 
2300 
2400 
2500 
2600 
2700 
2800 
2900 
3000 
3100 
320 ° 
3300 

3400 

3500 
3600 
3700 
3800 
3900 
4000 
4100 
4200 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
5000 
51 00 
5200 
5300 
5400 
5500 
5600 
6700 


FUNCTION DJMM(J,M1,M2*THETA) 

LOGICAL FIR ST/. TRUE./ , ODD 
REAL J,M1,M2 
DOUBLE PRECISION SUM 
COMMON /CGC0M/X( 200) 

DATA ZERO/5 .CE— 6/ 

DATA HAL FPl ,P I, T WOP I ,FOURP 1/1.5707963, 3.1415926, 6. 28 31 8 53, 12. 5 66371/ 
ODD(N)=MOD(N,2).FQ.l 

C THIS SUBROUTINE CALCULATES THF REDUCED ROTATION MATRIX <J , M i /R ( THE TA ) / J , M 2> 
IF (FIRST) Gl) TO 10 
C BEGIN CALCULATION 

1 D J MM=0 . 0 
N1M. l+( J + Ml) 

N2= 1 . 1 ♦ ( J-M2 ) 

N 3= 1 • 1 -M J— M 1 ) 

N4= 1 . 1 + ( J + M 2 ) 

IF <MIN0(N1,N2,N3,N4).LT. 1) RETURN 

P=X(N1 )+X(N2)+X(N3)+X( N4 ) 

FACTOR^ l .0 

C CHECK FOR STANDARD ANGLE 
ANGLF = THE T A 

IF ( ANGLE. LT. -ZERO) GO TO 7 
IF ( ANGLF.GT.PI ) GO TO 8 

2 IF ( ANGLE. LT.HALFPI ) GO TO 3 
AN GLF=Pl— ANGLE 

IF ( ODD ( N3 ) ) FACT OR =— FACT OR 
N 4=N 2 

3 N=N l +N3— 2 
N3=N1-N4 

C TEST FOR SPECIAL VALUES 

IF ( ANGLE. GT. ZFRO) GO TO 4 
D J MM=0 • 0 

IF (N3. FQ.O ) D JMM = F AC TOR 
RETURN 

C CONTINUF CALCULATION 

4 NMA X = MI NO ( N1 , N2 )-l 

NM IN = M AXO ( N 3,0 ) 

N P l =NM I N + 1 

IF ( ODD ( NM IN)) FACTOR=-FACTOR 
ANGLE=ANGLF/2.0 

IF (N.GT.O) F ACT 0R= FACT OR* (COS (ANGLE) ) **N 
T = T AN ( ANGLE ) 

N= 2*NM IN-N 3 
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5800 
5900 
6000 
6100 
6 200 
6300 
64-00 
6500 
6600 

6 7 00 
6800 
6900 
7000 
7100 

7 200 
7300 
7400 
7500 
7600 
7700 
7800 
7900 
R0C0 

00 

8 ? 00 
R300 
8400 
R 5 CO 

8 6 00 
870C 
R800 
8900 
9000 
91 00 
9200 
9300 
9400 


IF (N.GT.O) FACTOR=FACTOR*T**N 

T FRM= FACTOR 
SUM-T FRM 

IF (NMIN.FQ.NMAX) GO TO 6 
T=-T**2 

DO 5 N=NPl f NMAX 
ON E= ( N-N 1 ) *(N-N?) 

TW0=(N-N3)*N 
TFRM=TFRM*T*ONF 'TWO 
6 SIJ^SUM+TFRM 

6 N 1 =N l -NM f N 
N2=N2-NMIN 
N3 =N P 1 — N3 

Q=X(Nl>+X(N2)+X(N3)+X(N p L) 

0JMM=FX P<0. 5 * P-Q ) *SUM 
IF (ABSIDJMM) • 1 E* RO) DJMM=0.0 
RF TURN 

C PFOUCF ANGLF 

7 ANGIF=-ANGLE 

IF ( DDO ( N1-N4 ) ) FACTO R= -FACTOR 

8 ANGLF-AMOD(ANGLEtFOURPI) 

IF (ANGLE. LT.TWOPI ) GO TO 9 

ANGLE-ANGLF~TWO p I 

IF ( 000 ( N1 +N3 ) ) FACTOR=-FACTOR 

9 IF ( ANGLE. LT. PI ) GO TO 2 
ANGLE=TWO p I-ANGLE 

IF { ODD ( Nl + N2 ) ) FACTO R= -FACTOR 
GO TO 2 

C ARRAY X(N)=LOG(GAMMA(N) ) 


9500 


9600 

10 

FIRST=. FALSE. 

9700 


DO 11 N = 1 *2 00 

9800 


0 = N 

9900 

1 1 

X ( N)= ALGA^A ( Q ) 

1C00C 


GO TO 1 

10100 


END 



APPENDIX D 


SUBROUTINE FOR EVALUATING NINEJ COEFFICIENT 


100 
200 
300 
400 
500 
600 
700 
800 
900 
1000 
1100 
1200 
1300 
1400 
1500 
1600 
1700 
1800 
1900 
2000 
2100 
2200 
2300 
24 00 
2500 
2600 
2700 
2800 
2900 
3000 
3100 
3200 
3300 
3400 
3500 
3600 
3700 
3800 
3900 
4000 
4100 
4200 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
5000 
5100 
6200 
5300 
5400 
5500 


FUNCTION NINEJ<Ll,Sl,Jl,L2,S2»J2»L,StJ) 

LOGICAL ODD 

REAL LliSlt Jlf L2*S2t J2,LtSt J 

DOUBLE PRECISION SUM 

DATA ZER0fS0T4PI/5.0E-6»3* 5449077 / 

ODD(N)=M0DIN t 2).E0.1 

C THIS SUBROUTINE CALCULATES THE 9-J COEFFICIENT ( L 1 * SI * Jl » L2 , S2 » J2 t L « S t J ) 
C IT IS ASSUMED THAT ALL SELECTION RULES HAVE BEEN CHECKED PRIOR TO CALL 

P=AMIN1(L1+J*L2+S* J2+S1) 

Q= AMAX 1 (ABS <Ll~J)fABS< L2— S ) , ABS ( J 2-S 1 ) ) 

M=2.01*P 
N=2.01*Q 
M= M— N 

IF (ODD( M) .OR^M.LT.O) RETURN 
P=Q 

0 = 2* 0* P + 1.0 
N=M / 2+1 
SUM=0 • 0 
00 1 M=1 ♦ N 

SUM=SUM+Q*RACAH<Ll»L2fJ*StLfP)- 
1 *RACAH( J2.L2t Slf S, S2,P>- 

1 *R AC AH (Ll«SltJtJ?fJl«P) 

P=P + L.0 
1 0 = 0 + 2. 0 
N I NE J = SUM 
RETURN 

C THE FOLLOWING ENTRY IS FOR THE SPECIAL CASE Sl=S2=l/2 

C IT IS ASSUMED THAT ALL SELECTION RULES HAVE BEEN CHECKED PRIOR TO CALL 

ENTRY QLSJ( QLl tQJI tOL2tOJ2,QLtOS f QJ) 

M=ABS( OL-Q J ) +0 . 1 

N = Q S+0 • 1 

A=0 .0 

B= 0 • 0 

C=0.0 

0 = 0.0 

E= ( QL+Q J ) /2. 0 

IF (M.NE.O) GO TO 2 
A=RACAH(QLl f 0.5 t 0J tOJ2tQJl t QL2) 

IF (N.EQ.l) GO TO 2 
C=SQRT(4.0*QJ+2.0) 

GO TO 3 


2 D=4.0*QJ+2.0 
E = E + 0. 5 

B=D*RACAH(0L2*0.5 f 0Jl,0Jt0J2,E)*RACAH(QJl,0.5,0L2f0LtQLl,E> 
0=3.0*D*RACAH<0.5, 0.5, QL *QJ,1.0tE) 

3 NINEJ=( A-B)/ (C-0) 

RETURN 
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5600 
5 7 00 
5800 
5900 
6000 
6100 
6 ? 00 
6300 
6400 
6500 
6600 
6700 
6800 
6900 
7000 
71 00 
7200 
7300 
7400 
7500 
7600 
7700 
7800 
7900 
8000 
8100 
8200 
8300 
8400 
8500 
8600 
8700 


C THE FOLLOWING FNTRY IS FOR COMPUTATION OF 

C THF REDUCED MATRIX ELEMF NTS < LI * 1 /? , J1 / / TL S J // L 2 1 1 / 2 f J2> 

C IT IS ASSUMED THAT ALL SELECTION RULFS HAVE BEEN CHECKED PRIOR TO CALL 

FNTRY TLSJ(TL1,TJ1 ,TL2 ,TJ2*TL* TS,TJ ) 

NI NEJ=CG(TJ1 tO*5»TJt0«O»TJ2t0*5 I/SQT4PI 
N= T J + 0 « I 

IF (OnDCN)) N I NEJ=— N I NE J 
N=TS*(TJ-TL+2.1 ) 

IF CN.GT.O) GO TO 4 

C S=0 

N I N E J -N INEJ*SQRT( 2.0*TJ+1*0) 

RF TURN 

C S= 1 

4 A= ( T L+T J+l . 0 ) /2 * 0 

3= < TL 1-TJ l)*( 2.0*TJ1 + 1 .0) 

C= (TL2-TJ2 ) =M 2.0^TJ2+1 *0 ) 

GO TO (5*6*7) f N 

5 NINEJ=NINEJ*( B+C-A)/SQRT( A) 

RF TURN 

6 NI NEJ = NINEJ*( B-C } *SGR T { 2 # 0* A/TJ/(TJ+1*0) ) 

RETURN 

7 NI NEJ=NINEJ* ( 8+C+A)/S0RT ( A) 

RETURN 

END 
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